Information on the prioritizr package can be found at https://mran.microsoft.com/snapshot/2018-03-04/web/packages/prioritizr/vignettes/quick_start.html.
A really useful example for summed solutions is under “Selection Frequencies” here: https://cran.r-project.org/web/packages/prioritizr/vignettes/tasmania.html
For the solve() function, you must have a solver installed. For the the process with summing and multiple runs process, this must be gurobi. This requires following the instructions below.
install.packages('/Library/gurobi902/mac64/R/gurobi_9.0-2_R_3.6.1.tgz', repos=NULL) (if on a mac and downloaded the most recent version of gurobi, don’t download newest version of R! (4.0.0)).library(here)
library(janitor)
library(prioritizr) # marxan
library(sf) # spatial features
library(slam) # to use gurobi
library(gurobi) # solver
library(ggmap) # basemaps
library(patchwork)
library(tidyverse)
If you do not have any of these packages installed, run the following in the console: install.packages("package name")
Note: this will not work for installing gurobi. Instructions are in the introduction section.
All .csv files in this tutorial are unaltered from the .xlsx files provided in the R drive. They have only been renamed and saved as csvs. The “data” folder is identical to the MorroBay_data folder provided in the R drive, it has only been renamed.
spec <- read_csv("morro-bay-spec.csv") %>%
head(140) %>% # mine reads in extra blank rows - skip if yours does not
rename("amount" = "target") # target is called "prop" (relative) or "amount" (absolute)
pu <- read_csv("morro-bay-pu.csv") %>%
select(1:3) # mine reads in an extra blank column - skip if yours does not
puvsp <- read_csv("morro-bay-puvspr.csv") %>%
select(1:3) %>%
head(11849)
status <- read_csv("spec-name-status.csv") %>%
select(1:3) %>%
head(140)
parcels <- read_sf(dsn = here("data"), layer = "MorroBay_parcels") %>%
clean_names()
ggplot(data = parcels) +
geom_sf()
The function marxan_problem() in the prioritizr package gives us the “canned” approach that works for our purposes. If more customizations are desired, feel free to explore the problem() function instead.
marxan_1 <- marxan_problem(x = pu,
spec = spec,
puvspr = puvsp,
bound = NULL,
blm = 0)
marxan_1_problem <- marxan_1 %>%
add_gurobi_solver(gap = 0.15) %>%
add_pool_portfolio(method = 2, number_solutions = 100) # see method meaning under ?add_pool_portfolio
print(marxan_1_problem)
marxan_1_soln <- solve(marxan_1_problem)
## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 140 rows, 670 columns and 11849 nonzeros
## Model fingerprint: 0xd2c5bccb
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
## Matrix range [1e+00, 1e+00]
## Objective range [2e+00, 5e+06]
## Bounds range [1e+00, 1e+00]
## RHS range [3e-01, 2e+02]
## Found heuristic solution: objective 4.238988e+07
## Presolve removed 116 rows and 205 columns
## Presolve time: 0.01s
## Presolved: 24 rows, 465 columns, 3066 nonzeros
## Variable types: 0 continuous, 465 integer (465 binary)
## Presolve removed 1 rows and 0 columns
## Presolved: 23 rows, 465 columns, 3058 nonzeros
##
##
## Root relaxation: objective 2.840941e+07, 14 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 2.840941e+07 2.8409e+07 0.00% - 0s
##
## Optimal solution found at node 0 - now completing solution pool...
##
## Nodes | Current Node | Pool Obj. Bounds | Work
## | | Worst |
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 - 0 - 2.8409e+07 - - 0s
## 0 0 - 0 - 2.8409e+07 - - 0s
## 0 2 - 0 - 2.8409e+07 - - 0s
##
## Explored 1958 nodes (383 simplex iterations) in 0.13 seconds
## Thread count was 1 (of 4 available processors)
##
## Solution count 100: 2.84094e+07 2.84097e+07 2.842e+07 ... 3.34128e+07
##
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.840941400000e+07, best bound 2.840941400000e+07, gap 0.0000%
marxan_1_ssoln <- marxan_1_soln %>%
mutate(sum = rowSums(.[6:105])) %>%
select(id, cost, status, locked_in, locked_out, sum)
hist(marxan_1_ssoln$sum)
parcels_marxan_1 <- inner_join(parcels, marxan_1_ssoln, by = "id")
ggplot(data = parcels_marxan_1) +
geom_sf(aes(fill = sum),
color = "white",
size = 0.05) +
scale_fill_binned(low = "slategray2",
high = "navy") +
labs(fill = "Summed \nSolution") +
theme_minimal()
This step is optional, but will make your map look better and is good to learn for further spatial analyses in R. Unfortunately, the options for basemaps in R are limited. Here, I have used the ggmaps package because it is the most accessible, however, it is still not very accessible compared to other R packages.
Using ggmap requires getting an API key, which is a bit of a hassle, but worth it. Instructions are here: https://cran.r-project.org/web/packages/ggmap/readme/README.html
For a ggmap cheatsheet, including the different basemaps in ggmaps, see: https://www.nceas.ucsb.edu/sites/default/files/2020-04/ggmapCheatsheet.pdf
morrobay <- get_map(location = c(lon = -120.7665, lat = 35.335),
zoom = 12,
maptype = "terrain-background", # - background means no references, omit if want references
source = "google")
all_spec <-
ggmap(morrobay) +
geom_sf(data = parcels_marxan_1,
aes(fill = sum),
color = "white",
size = 0.1,
alpha = 0.85,
inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326)) +
scale_fill_gradient(low = "tan",
high = "red3") +
labs(title = "All Species",
fill = "Summed \nSolution",
x = NULL,
y = NULL) +
theme_minimal()
all_spec
end_status <- status %>%
mutate("endangered" = case_when(
str_detect(status, pattern = "endangered") == TRUE ~ "yes",
str_detect(status, pattern = "threatened") == TRUE ~ "yes",
T ~ "no")) %>%
filter(endangered == "yes")
end_spec <- merge(end_status, spec, by = "id") %>%
select(id, amount, spf, name.x) %>%
rename("name" = "name.x")
puvsp_id <- puvsp %>%
rename("id" = "species")
end_puvsp <- merge(end_status, puvsp_id, by = "id") %>%
rename("species" = "id")
marxan_end <- marxan_problem(x = pu,
spec = end_spec,
puvspr = end_puvsp,
bound = NULL,
blm = 0)
marxan_end_problem <- marxan_end %>%
add_gurobi_solver(gap = 0.15) %>%
add_pool_portfolio(method = 2, number_solutions = 100)
print(marxan_end_problem)
marxan_end_soln <- solve(marxan_end_problem)
## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 18 rows, 670 columns and 970 nonzeros
## Model fingerprint: 0x6e88f508
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
## Matrix range [1e+00, 1e+00]
## Objective range [2e+00, 5e+06]
## Bounds range [1e+00, 1e+00]
## RHS range [9e-01, 1e+02]
## Found heuristic solution: objective 4.174920e+07
## Presolve removed 14 rows and 203 columns
## Presolve time: 0.00s
## Presolved: 4 rows, 467 columns, 495 nonzeros
## Variable types: 0 continuous, 467 integer (467 binary)
## Found heuristic solution: objective 3.534249e+07
## Presolved: 4 rows, 467 columns, 495 nonzeros
##
##
## Root relaxation: objective 2.677119e+07, 1 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 2.677119e+07 2.6771e+07 0.00% - 0s
##
## Optimal solution found at node 0 - now completing solution pool...
##
## Nodes | Current Node | Pool Obj. Bounds | Work
## | | Worst |
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 - 0 - 2.6771e+07 - - 0s
## 0 0 - 0 - 2.6771e+07 - - 0s
## 0 2 - 0 - 2.6771e+07 - - 0s
##
## Explored 4914 nodes (1716 simplex iterations) in 0.24 seconds
## Thread count was 1 (of 4 available processors)
##
## Solution count 100: 2.67712e+07 2.67732e+07 2.67736e+07 ... 3.14836e+07
##
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.677119300000e+07, best bound 2.677119300000e+07, gap 0.0000%
marxan_end_ssoln <- marxan_end_soln %>%
mutate(sum = rowSums(.[6:105])) %>%
select(id, cost, status, locked_in, locked_out, sum)
hist(marxan_end_ssoln$sum)
parcels_marxan_end <- inner_join(parcels, marxan_end_ssoln, by = "id")
ggplot(data = parcels_marxan_end) +
geom_sf(aes(fill = sum),
color = "white",
size = 0.05) +
scale_fill_binned(low = "slategray2",
high = "navy") +
labs(fill = "Summed \nSolution") +
theme_minimal()
end_spec <-
ggmap(morrobay) +
geom_sf(data = parcels_marxan_end,
aes(fill = sum),
color = "white",
size = 0.1,
alpha = 0.85,
inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326)) +
scale_fill_gradient(low = "tan",
high = "red3") +
labs(title = "Endangered & Threatened Species",
fill = "Summed \nSolution",
x = NULL,
y = NULL) +
theme_minimal()
end_spec
This step isn’t necessary, but is great if you want to present your maps together. Patchwork uses PEMDAS to combine ggplots together into one graphic. For example, plot_1 + plot_2 will make an image of the plots side by side, whereas plot_1 / plot_2 will make an image of plot_1 over plot_2. For more information, see https://github.com/thomasp85/patchwork
# Make the all_spec graph without a legend so that the combined image has only one legend
all_no_lgnd <- all_spec +
theme(legend.position = "none")
# Combine graphs
spec_graphs <- all_no_lgnd + end_spec
spec_graphs
# Save graphs as .png to add to your report
ggsave("spec-graphs.png")